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ABSTRACT.  This  paper  presents  the  survey  of  the  basic  ideas  and  results  of 


the  h,  p  and  h-p  versions  of  the  finite  element  method. 
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1 .  INTRODUCTION 

Finite  element  method  is  today  the  major  tool  in  computational 
mechanics.  Many  research  and  commercial  programs  are  used  in  practice.  More 
than  40  thousand  papers  of  theoretical  and  computational  characters  have  been 
published. 

The  software  of  the  finite  element  method  is  complex  and  is  influenced 
by  the  theory,  implementational  aspects,  computational  experience  and 
engineering  needs.  It  is  influenced  by  the  hardware  development  and  the 
relation  between  the  computer  and  manpower  cost  as  well  as  the  class  of 
problems  to  be  solved  by  the  engineering  community.  The  finite  element  is  a 
tool  for  obtaining  the  data  of  interest  with  the  accuracy  in  a  prescribed 
range  in  an  effective  way  (taking  into  account  the  computer  and  manpower 
cost. )  The  existence  of  a  reliable  and  accurate  quantitative  a-posteriori 
error  estimation  for  anv  data  of  Interest  should  be  one  of  the  major  aspects 
of  assessing  the  finite  element  method  in  general  and  a  code  in  particular. 

The  various  adaptive  approaches  together  with  users  interaction  are  essential 
for  the  effectiveness  of  the  approach. 

The  majority  of  the  results  and  finite  element  programs  in  structural 
mechanics  are  based  on  the  classical  h-version.  Nevertheless,  during  the  last 
ten  years  the  p-version  and  the  h-p  version  was  developed,  major  theoretical 
results  achieved  and  few  large  scale  program  developed.  Especially,  we 
mention  here  the  program  STRIPE  (Aeronautical  Research  Institute  of  Sweden), 
Applied  structure  (Rasna  Corp. ,  CA,  USA),  PHLEX  (Computational  Mechanics,  TX, 
USA)  and  MSC/PROBE  MacNeal  Schwendler,  CA,  USA).  (We  note  that  less  than 
one  hundred  papers  addressing  the  p  anH  h-p  versions  have  been  published. ) 

Comparison  of  the  methods  and  codes  based  on  the  h,  p  and  h-p  versions 
is  a  complex  task  and  depends  on  the  criteria  used.  It  seems  that  one  of  the 
ma lor  criterion  is  the  existence  of  a  reliable  quantitative  error  estimation 
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of  any  data  of  interest — not  only  for  example  measured  by  the  energy  norm.  In 
addition,  the  adaptive  approaches  are  imperative,  together  with  the  robustness 
of  the  method  and  its  effectiveness  in  the  sense  of  computer  and  manpower 
cost.  Theory  of  the  method  has  to  be  well  developed,  so  that  it  reliably 
serves  to  the  understanding  of  the  method  and  its  features. 

In  this  paper  we  will  focus  on  very  basic  aspects  of  the  h,  p  and  h-p 
versions  of  the  finite  element  method,  especially  with  respect  to  the 
theoretical  understanding.  Me  will  show  that  it  is  possible  to  obtain  very 
useful  understanding  from  the  detailed  mathematical  analyses  of  a  simple  one 
dimensional  problem.  The  mathematical  theorems  properly  used  are  very 
important  guidelines  for  the  software  in  2  and  3  dimensions  too.  We  will 
present  also  basic  results  in  2  dimensions  and  highlight  the  similarity  and 
differences  between  the  one  and  two  dimensional  cases.  In  this  paper  we  will 
focus  on  aspects  and  examples  of  most  simple  nature,  but  still  keeping  certain 
essentials  of  the  general  cases  intact.  We  will  emphasize  the  aspects  of  the 
E  and  h-p  versions  of  the  finite  element  method.  For  the  survey  of  the 
practical  aspect  of  the  p-version  in  the  context  of  solving  complex  3 
dimensional  program  based  on  program  STRIPE,  we  refer  to  [2]  [3],  For  the 
survey  of  the  basic  theoretical  results,  we  refer  to  [20]. 

In  Section  2  we  will  address  the  h,  p  and  h-p  versions  in  one  dimension, 
present  basic  very  detailed  theoretical  results  and  make  comments  of  their 
meaning  and  importance  in  more  general  setting  of  the  adaptivity,  etc. ,  with 
analogs  in  2  and  3  dimensions.  Section  3  will  elaborate  on  the  two 
dimensional  problems. 

In  this  paper  we  will  restrict  ourselves,  for  simplicity,  only  to  the 
energy  norm  accuracy  measure,  although  it  has  relatively  very  limited 
practical  engineering  importance  in  the  sense  we  defined  above,  but  I  will 
still  allow  us  to  show  many  basic  features  of  the  general  case  (discusses  in 
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[2]  [3]).  Numerical  examples  will  illustrate  the  use  of  the  presented 
theoretical  result  in  practicai  context. 

2.  THE  FINITE  ELEMENT  METHOD  IN  ONE  DIMENSION 
2.1.  Formulation  of  the  problem 

As  said  in  the  introduction,  the  one  dimensional  problem  is  a  good  way 
to  explain  the  essentials  of  the  h,  p  and  h-p  versions. 

We  will  consider  the  following  simple  one  dimensional  model  problem 

(2.1.1a)  -u"  =  f  in  I  =  (0,1) 

(2.1.1b)  u(0)  =  u(l)  =  0 

with  the  exact  solution 

(2.1.2)  Uq(x)  =  (x“-x).  a  >  1/2. 

Let  H^(I)  be  the  standard  Sobolev  space  and 

hJ(I)  =  (u  e  hJ(I)|u(0)  =  u(l)  =  0>. 

Further,  let 

(2.1.3)  B(u,v)=ru'v'dx 

11  1  /O 

be  the  bilinear  form  on  H^d)  x  H^d)  and  let  =  (B(u,u)) 

Obviously,  ||•||p  is  equivalent  on  H^d)  with  the  Sobolev  norm  1*1  .  We 

H  (I) 

will  assume  that  a  >  1/2  to  get  ||’Jq||£  *• 

The  solution  (2.1.2)  is  a  very  good  model  for  the  two  dimensional 

problem  where  the  domain  has  corners  (e.g.,  cracked  domain,  etc.)  or 

interfaces  and  of  3  dimensional  problems  too. 

B  1 

In  two  dimensions  u  =  r^v(0)  €  H  (£1)  for  B  >  0.  This  and  other 
reasons  show  that  the  results  in  two  dimensions  with  the  singularity  of  the 
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type  r^(p(6)  should  be  compared  with  the  one  dimensional  case  when  a  = 

^  +  1/2. 

Let 

A  ;=  {0  =  <  xf  <  •••  <  =  1} 

0  1  m(A) 

be  the  mesh  on  I  and  =  (xf  .,x^),  i  =  l,...,m(A),  h^  =  I  1^1  =  xf-x^  , 

h(A)  =  max  h^,  i  =  1 ,2, . . . ,m(A) .  The  points  x^  will  be  called  the  nodal 

points  and  the  elements.  Let  further  p  =  |p^ . ^mCA)]’  ^i  ~ 

integer  be  the  element  degree  vector.  We  will  use  the  notation  E  =  (A,p) 
which  will  characte'ize  the  basic  finite  element  method.  Let 

S  =  S(Z)  =  {u  €  Hq(I)  u[j  is  a  polynomial  of  degree  p^> 

i 

We  will  call  N(Z)  =  dim  S(Z)  the  number  of  degrees  of  freedom. 
Finite  element  solution  u^  e  S(Z)  is  then  defined  so  that 

(2.1.4)  B(u2.,v)  =  B(Uq,v),  V  V  e  S(Z). 


We  will  denote  the  error  by  e  =  e^.  =  u^.  -  u  and  will  be  interested  in 
the  accuracy  measure  ||*||£ 

Let  us  now  define  the  h,  p  and  h-p  versions  of  the  finite  element 
method.  To  this  end,  assume  that  a  sequence  Z^^  with  N(Zj^)  =  oo 

be  given  (resp.  constructed  by  an  adaptive  approach).  Then  computation  of 
u_  will  be  called 

'  fA.  A.. 

a)  The  h-version:  if  Zj  =  (Aj.p^),  where  p^  =  Ipg'’ . ^0^1’ 

i.e.,the  degree  of  the  elements  are  fixed  and  the  mesh  is  changed 
(refined).  This  is  the  classical  version  of  the  finite  element 
method. 

b)  The  D-version:  if  Z.  =  (A.d.).  where  =  ID^  . ,.,) 

-  J  j  J  t  J.l  ^j,m(A) 
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either  (i.e.,  uniform  p)  or  p^  .  5*  p^  „  if  i*  i 

J.  1  J  J, 1  J.« 

(i.e.,  nonuniform  degrees),  i.e.,  mesh  is  fixed  and  degrees  are 
increased. 

d)  The  h-D  version:  if  Z.  =  (A..d.)  where  both,  the  meshes  and  the 
- -  J  J  J 

degrees,  are  changed. 

We  will  associate  to  the  solution  u^  a  computable  a-posteriori  error 


estimator 


|e^||^.  Usually  the  error  estimator  is  computed  by  elemental  error 


indicator  7)(I.)  with 
-  1 


r(A)  -  .  ,1/2 

i=l 


We  will  assume  that  the  goal  of  computation  is  to  design  Z  so  that  ~ 

lezlg  -  t||u||^  where  t  is  an  a-priori  given  tolerance  (say  1%  or  5%).  We 
will  not  address  here  in  detail  the  concrete  form  of  the  estimation. 

2.2.  The  h-verslon  of  the  finite  element 

Obviously  we  have  for  the  h-version  N(Z)  =  mp  -  1.  The  first 
theoretical  question  is  what  is  the  best  accuracy  the  h-versior  could  provide 
among  all  meshes.  We  have 

Theorem  2.2.1  (36].  Let  a  >  1/2  be  a  non-integer,  and  u^  be  given 
by  (2.1.2),  then  there  is  a  constant  C  =  C(a,p)  >  0  such  that  for  any  mesh 

-P 

(2.2.1)  |]ej.||^  £  CN  □ 

This  theorem  shows  that  (asymptotically)  it  is  not  possible  to  expect 
anything  better  than  the  rate  N  ^  when  elements  of  degree  p  is  used.  This 
does  not  mean  that  this  accuracy  will  be  achieved  for  any  mesh. 

The  question  is  whether  for  a  proper  sequence  of  meshes  we  get 
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I^zIe  ^  ■ 


Hence,  the  problem  arises,  what  is  the  best  theoretical  mesh.  To  this  end  we 
introduce  the  notion  of  the  grading  function.  Function  g(x),  0  £  x  s  1  is 
called  a  mesh  grading  function  if  the  nodal  points  of  the  mesh  are  such  that 


^i  =  S 


(i)-  ‘  . "  = 


m(A) 


The  mesh  will  be  called  radical  if  g(x)  =  x^.  We  shall  assume  that  the 
grading  function  g(x)  satisfies  the  following  conditions: 

G(l):  g(0)  =  0,  g(l)  =  1, 

G(2);  g  is  continuous  and  strictly  increasing, 

G(3):  g  €  C^(0,1)  n  C°[0,1]. 

Then  we  have 

Theorem  2.2.2  [36].  Among  all  grading  functions  g(t)  satisfying 
Assumption  G(l)  -  G(3) 

®opt  ^  a-1/2 

is  the  optimal  one.  Precisely  with  this  grading  function 


(2.2.2) 


pA 


lim  mP|ej.||g.  =  C(a,p)  f  — ?  1 

m-^  '  ' 


where 

(2.2.3) 


,  ar(a)  sin  nal  r(p-a+l) 
Cia,b)  =  - - L  - c 


Vn  4^v'2p+l  r(p+^) 


attains  the  minimum.  We  see  that  the  optimal  mesh  is  radical  one  and  gives 
maximal  possible  rate  N 


Theorem  2.2.3  [36] : 
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then 


a) 


If 


(2.2.4) 


lim 

m-)oo 


'E"E 


1 

C(a,p)3 


✓(2a-l)p-2p 


b)  If  p 

cr 


(2.2.5) 


then 


lim 

m-xn 


Vln  m 


1 

C(a,p)/3 


c) 


If 


then 


(2.2.6) 


lim  m^||ej,|^  =  C^(a,^,p) 

m-xn 


where  C(a,p)  is  defined  by  (2.2.3)  and  0  <  Cj(a,/3,p)  <  »  has  a  more 
complicated  expression.  □ 

We  remark  that  0=1  corresponds  to  the  uniform  mesh.  Theorems 
2.2. 1-2. 2. 3  Indicate 

a)  The  maximal  achievable  rate  of  convergence  is  N  ^  and  this  rate 
can  be  achieved  by  a  proper  mesh  also  for  nonsmooth  solution  of  the 
type  (2.2.1).  This  optimal  mesh  depends  on  the  strength  of  the 
singularity  and  the  used  degrees  of  the  elements.  Increasing  the 
degrees  of  elements  leads  to  the  strengthening  of  the  refinement. 

b)  The  optimal  mesh  is  the  radical  one.  Underref inement  is  more 
"dangerous"  than  overref inement,  because  it  can  significantly  slow 
down  the  convergence. 

We  have  addressed  the  method  and  the  optimal  mesh  only  in  the  case  when 
the  solution  is  given  by  (2.2.1).  In  general,  the  mesh  has  to  be  constructed 
by  an  adaptive  procedure.  Two  types  of  adaptive  mesh  construction  are  used 
(in  one  or  2,  3  dimensions). 

a)  A  sequence  of  meshes  in  constructed  by  consecutive  refinement  or 
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derefinement  of  the  particular  elements.  The  goal  for  refinement 
or  derefinement  is  to  make  the  (elemental)  error  indicators  tj 
approximately  equal.  For  a  detailed  theoretical  analysis  of  this 
approach  in  one  dimension,  we  refer  to  [23].  The  final  mesh  is 
accepted  when  the  error  estimator  indicates  the  acceptable  accuracy. 

b)  The  initial  usually  crude  mesh  is  used  for  the  computation  of  the 

basic  characteristic  of  the  solution.  Using  these  characteristics  a 
new  mesh  is  constructed  with  the  aim  to  get  acceptable  accuracy.  If 
the  accuracy  is  not  achieved,  the  process  is  repeated.  (See.  e.g. 

[6] [39] [40] [57] [591 ) 

We  will  call  the  adaptive  approach  based  on  a)  the  local  approach  and  if 
it  is  based  on  b)  the  global  approach. 

The  mentioned  theoretical  results  give  an  important  insight  into  the 
design  of  the  adaptive  procedure  and  judging  its  performance  when  used  on  the 
benchmark  problem  with  the  solution  (2.2.1).  The  notion  of  the  grading 
function  can  be  used  in  general  and  especially  in  the  global  adaptive 
approach.  See  [39]  [40]. 

The  mesh  construction  is  very  influenced  by  the  implementational 
aspects.  For  example,  in  the  local  adaptive  approach  often  only  the  halving 
of  particular  elements  is  used.  Accuracy  and  reliability  control  is  a  very 
important  task  in  the  practice  of  finite  element  method.  Various  mathematical 
and  heuristic  methods  for  1,  2  and  3  dimensional  problems  were  suggested.  We 
refer  here  as  examples  to  [ 1 ] [5] [6] [ 14] [ 16] [ 17] [ 18] [ 19] [24] [25] [30] [31 ] [32] 
[34] [40] [41] [47] [48] [53] [54] [55] [56] [57] [58] [60] [51]  and  the  survey  [17]. 
Various  adaptive  codes  based  on  the  h-version  were  written  (see.  e.g.  [26] 

[33]  [44]  [51]  [59]).  The  error  estimates  are  in  the  h-version  usually  for 
the  energy  norm  error.  The  a-posteriori  estimates  for  other  data  as  stresses 
in  a  point,  the  stress  intensity  factor,  etc. ,  are  almost  impossible  to  get  by 
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the  approaches  used  for  the  energy  norm  measure. 


2.3.  The  p-version 

In  the  p-version  the  mesh  A  is  fixed,  i.e.,  the  mesh  is  made  usually 
by  the  user  (based  on  his  experience  and  various  basic  rules)  and  the  accuracy 
is  obtained  by  (adaptively)  increasing  the  degrees  of  elements.  The  degrees 
could  be  uniform  or  nonuniform.  The  error  estimator  S  is  usually  based  on 
comparison  of  computed  data  for  the  sequence  of  used  degrees  and  an 
extrapolation  procedure.  See,  e.g.  [54,55,56]  for  details. 

In  2  and  3  dimensions,  the  user  identifies  the  singular  areas  as  corners 
and  edges,  and,  for  example,  constructs  locally  refined  meshes  in  these 
singular  areas,  and  then,  by  a  general  mesh  generator,  connects  (complements) 
these  local  meshes  to  the  mesh  of  the  entire  domain  or  constructs  first  the 
mesh  and  then  refines  it  in  critical  areas.  The  degree  of  the  elements 
which  leads  to  the  desired  accuracy  is  then  adaptively  determined.  The 
singularity  refinement  is  made  by  the  "layers"  of  the  elements  around  the 
singularity.  (See  examples  of  two  dimensional  meshes  of  this  type  in  Section 
3  (Fig.  3.2.1)).  A  minimum  of  2-3  layers  is  practically  necessary  because  for 
accurate  stress  computation  in  a  point  at  least  2  layers  have  to  separate  the 
point  from  the  singularity  point.  The  p-version  allows  to  use  relatively  very 
small  number  of  elements. 


Let  us  now  present  typical  results  showing  the  features  of  the  p-version 


when  applied  to  the  solution  (2.1.2).  To  this  end  let 


(2.3.1) 


rj(I^)  =  inf  ||u  -  u|l 

u=P  (I.)  ^  Ed.) 

pi  1 


Z"  A 
^  (d) 
1 


where  P  (1^^)  is  the  class  of  all  polynomials  of  degree  p.  We  mention  that 


(2.3.1)  is  in  our  case  the  error  of  the  finite  element  solution  and 


1=1 


9 


Theorem  3.3.1  shows  that  in  the  first  element  the  rate  is  only 
algebraic,  while  in  all  others  it  is  exponential.  It  shows  also  that  the 
accuracy  in  the  first  element  is  essentially  governed  by  its  size  (h^  = 

b)  while  the  accuracy  of  the  other  elements  is  governed  by  their  degree 
(because  the  exponential  rate).  We  see  that  the  error  ij(I^),  i  a  2  will  be 
essentially  equal  (for  uniform  p)  in  all  the  elements  when  the  mesh  is 

geometric  with  the  factor  q,  i.e.,  x.  ,  =  qx^,  x^,,,  =  1.  Hence,  a 
geometric  mesh  is  a  good  choice  (see  Section  2.2.4  for  additional  arguments). 

The  question  how  many  elements  and  which  factor  q  should  be  selected 
is  very  important  for  practical  computation.  Small  q  and  large  m(A)  can 
lead  to  the  "overrefined"  mesh  for  p  small  because  the  error  ||®||£  is 
governed  by  the  error  in  the  largest  element.  For  high  p  the  small  factor 
q  is  preferable  because  the  error  is  essentially  governed  by  the  elements 
close  to  the  singularity.  The  overref inement  is  advisable  because  when  small 
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p  (which  is  adaptively  determined)  leads  to  the  desired  accuracy,  the 
computation  is  cheap  and  overref inement  does  not  essentially  increase  the 
cost.  On  the  other  hand,  if  the  desired  accuracy  needs  high  p,  the 
underref inement  could  be  very  costly.  Hence,  the  user  designing  a  mesh  has  to 
take  a  "cost  risk"  although  he  will  get  always  the  required  accuracy  by 
adaptive  solution  of  the  degrees. 

To  illustrate  these  points,  consider  the  mesh  with  2  elements  = 
(0,p),  I„  =  (p,l)  and  compute  (see  (2.3.2)) 

1 

“'2  1 

(2.3.5)  T,^(I^)  =  p6  ^ 

P 


(2.3.6) 


Tj.dz) 


2 


P 

for  various  p,  a  and  5.  Note  that  in  (2.3.5)  and  (2.3.6)  we  did  not  include 

the  constants  ^^(a)  and  C^(a)  which  have  the  factor  |sin  7ra| .  To  model 

1/2 

the  case,  which  we  will  study  in  Section  3  (u  »  r  )  we  will  use  a  = 

1,2,3.  This  has  to  be  understood  as  the  limiting  case  a  -  l-e,2-e,3-e  with 

*  ct 

G  ^  0  or  when  the  approximate  solution  is  u  =  x  Ig  x. 

2  2 

In  Table  2.3.1,  we  report  t7,(I^)  and  as  function  of  p,  p  and 

a  =  1,2,3.  We  see,  for  example,  that  for  a  =  1  the  error  is  governed  by  the 
error  in  for  p  a  2,  and  hence  p  is  "too  large",  i.e.,  the  mesh  is 

underrefined.  F'or  a  =  2,  the  mesh  is  "overrefined"  when  p  =  0. 1  and  p  £ 
8,  and  for  0.2  when  p  s  5.  Noting  the  magnitude  of  the  errors,  we  see 
that  the  overref inement  makes  no  harm  because  the  accuracy  is  achieved  for 
small  p,  and  computation  is  cheap.  From  (2.3.5)  and  (2.3.6)  as  well  as  from 
Table  2.3.1,  we  see  that (for  the  fixed  mesh)  for  smaller  p  the  rate  is 
exponential  and  for  p  larger  when  is  error  governing  interval,  the  rate 
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TABLE  3.1.  The  errors 


a  =  1 


=  .  15 


=  0.5 


2.42-1 

1.50-1 

1.65-1 

5.01-1 

1 . 63-2 

3.75-2 

8.08-3 

1.25-1 

1 . 96-3 

1 . 66-2 

7.00-4 

5.55-2 

2.98-4 

9.37-3 

7.69-5 

3. 12-2 

5. 15-5 

6.00-3 

9.59-6 

2.00-2 

9.65-6 

4. 17-3 

1 . 30-6 

1.38-2 

1.91-6 

3.06-3 

1 . 86-7 

1 . 02-2 

3.95-7 

2.34-3 

2.78-8 

7.81-3 

a  -  2 


=  0.2 


=  0.5 


B 

1 . 00-3 

1.96-1 

3.83-3 

1.19-1 

1.25-1 

2 

1 . 56-5 

3.31-3 

5.27-5 

1 . 46-3 

1 . 95-3 

3 

1 . 37-6 

1.77-4 

4.63-6 

5.62-5 

1.71-4 

B 

2.44-7 

1.51-5 

8.24-7 

3.47-6 

3.05-5 

5 

6.40-8 

1 . 66-6 

2. 16-7 

2.77-7 

8.00-6 

6 

2. 14-8 

2. 17-7 

7.23-8 

7.23-8 

2.68-6 

B 

8.50-9 

3.65-8 

2.81-8 

2.74-9 

1 . 06-6 

8 

3.81-9 

5.01-9 

1 . 29-8 

3. 14-10 

4.76-7 

12 


TABLE  3.1.  The  errors  { I  ). 

(Continued)  ^  ^ 


I 

I 


is  algebraic.  The  transition  point  between  the  exponential  and  algebraic  rate 
is  when  the  errors  in  and  are  equal.  We  will  see  the  same 

features  in  2  dimensions. 

Theorem  2.3.1  and  Table  2.3.1  show  also  that  nonuniform  degree 
distribution  is  optimal.  If  the  mesh  is  underrefined,  then  high  p  in 
and  low  in  is  optimal  and  for  overrefined  mesh  the  opposite  holds.  The 

adaptive  code  for  nonuniform  degrees  constructs  such  distribution.  Once  more 
we  see  that  overref inement  is  preferable.  The  value  of  p  ~  0.15  is,  as  we 
will  see  in  the  next  section,  optimal  for  all  a  when  nonuniform  p  is  used 
and  practically  is  good  also  for  p  uniform.  The  same  effects  hold  in  2 
dimensions  as  will  be  seen  in  Section  3. 

As  was  said  above,  the  mesh  design  is  very  important.  We  refer  to  [2] 

[3]  [54]  [55]  [56]  for  practical  rules  for  a  design  as  a  "good"  mesh.  Let  uj 


1.02-13  5.41-8 


1.65-13  4.88-9 


3.59-14  5.23-10 


9.33-15  6.33-11 


1 . 28-9 


7.24-11 


7.77-12 


1.26-12 


2.69-13 


7.07-14 


3.54-12 


1.10-10  5.08-18 


2.91-11  6.72-20 
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underline — as  we  have  seen — that  meshes  for  the  p-version  have  different 
character  than  that  of  the  h-version  with  fixed  degree  p.  Nevertheless, 
sometimes  an  adaptive  construction  of  the  mesh  for  the  h-version  with,  say  p 
=  3,  could  be  possibly  used  for  the  p-version  although  it  is  far  from  the 
optimal  [59].  This  could  be  seen  from  comparison  of  the  formulae  for  the 
optimal  mesh  for  the  h-version  and  p  *  3  discussed  in  Section  2.2  and  the 
formulae  of  Section  2.3.  The  error  estimator  based  on  the  extrapolation  is 
very  accurate  and  robust.  We  will  see  it  more  in  2  dimensional  example. 

The  first  theoretical  results  for  the  p-version  in  2  dimensions  are  in 
[22] ,  for  the  three  dimensional  problems  in  [28]  [29] .  For  the  survey  of  the 
results,  we  refer  to  [20].  The  p-version  is  related  to  the  spectral  element 
method.  See,  e.g.  [43]. 

2.4.  The  h-p  version 

In  this  section  we  will  discuss  theoretical  questions  about  the  h- 
version.  The  theoretical  results  will  also  give  an  additional  insight 
into  the  desirable  mesh  design  for  the  p-version. 

The  first  question  is  about  the  lower  bound  of  the  error. 

Theorem  2.4.1  [36].  For  any  mesh  A  with  m(A)  elements  and  any 
p  =  j  degree  vector  we  have 

'In 

(2.4.1)  ||e_|l  2:  C(a)  — - 

II  2^  H  /■— 

where 

(2.4.2)  p  =  a  -  1/2,  q^  =  {.VZ  -  1)^  «  0.17. 

(Here  once  more  the  solution  i^  is  given  in  (2.2.1)).  d 

The  next  question  is  whether  the  error  (2.4.1)  can  be  achieved. 

Theorem  2.4.2  [36]:  Let  the  mesh  A  be  geometric  with  the  factor  q. 
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„  A  m(A)-i  .  jii. 

i.e.  Xq  =  0,  Xj^  =  q  ,  1  =  1 ,2,  .  .  .  ,m(AJ  and  let  = 

(where  (•]  means  the  integral  part).  Then  we  have 

1.  if  s  >  Sq,  then 


[ l+s( i-1 ) ] 


(2.4.3) 


|ej.||  »  C(a,q,s)q 


(a-A)v^N7s 


2.  if  s  <  Sq,  then 


(2.4.4) 


3.  if  s  =  Sq,  then 


(2.4.5) 


(2.4.6) 


->/(a-^)N  'JZ  £n  q  £n  r 


and 


|ej.||  a  C(a,q)e 


r  =  -  (a  -  1/2)  ^  . 


Furthermore,  the  optimal  geometric  mesh  and  linear  degree  vector  combination 


is  given  by 
(2.4.7) 


%pt  = 


s  .  =  2a  -  1 

opt 


In  this  case 


(2.4.8) 


'Z"E 


CU)l{V2  -  1)^] 


)N 


In  (2.4.3)  and  (2.4.7)  ~  means  equivalency  with  the  constants  depending  on 

the  (a,q,s),  (a,q)  and  a,  respectively.  d 

Comparing  (2.4.1)  and  (2.4.3)  we  see  that  up  to  the  algebraic  term 
vS'  in  (2.4.1)  the  error  for  the  geometric  mesh  is  optimal  (i.e.,  the 
geometric  mesh  can  be  in  this  sense  understood  to  be  optimal). 

Theorem  2.4.1  and  (2.4.2)  show  that  geometric  mesh  with  q  «  0.17  is 
optimal  independently  of  the  strength  of  the  singularity  a  provided  that 
degrees  of  the  elements  are  not  constants  and  their  distribution  depend  on  the 
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strength  of  the  singularity.  For  the  uniform  degree  distribution  and 
geometric  mesh,  we  have 


Theorem  2.4.3  [36]  Let  A  be  a  geometric  mesh  with  the  factor  q  and 
m  elements  and  uniform  degrees  p  =  sm(A).  Then  we  have 


1)  If  s  >  s. 


(2.4.9) 


~  C(a.q) 


(a-^) vNTs 


2)  If  s  <  s. 


(2.4. 10) 


3)  If  s  =  Sq,  then 


(2.4.11) 


i®ziE  C(a.q) 


v(a~j  )  N^n  r  £n  q 


is  optimal  for  given  q  with 


(a-2)^n  q 


<r  =  min(2a-l,a). 


The  optimal  combination  is  given  by 


(2.4. 12) 


“I  '  %pt  = 


(2.4. 13) 


s  =  s  .=  2a  -  1 
opt 


Theorem  2.4.3  shows  that  the  factor  q^  of  the  geometric  mesh  is  optimal  for 
uniform  mesh  but  the  exponential  rate  of  the  convergence  is  not  optimal.  The 
exponent  is  by  the  factor  smaller.  We  see  also  that  the  increase  of  the 
degrees  in  analogous  as  in  the  case  of  nonuniform  degrees. 


Theorem  2.4.3  does  not  apply  that  for  the  uniform  degree  the  geometric 
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mesh  is  optimal.  In  fact,  we  have  ween  in  Section  2.2  that  for  every  degree 
p  the  radical  mesh  is  optimal.  Hence,  we  can  ask  the  question,  what  rate  of 
convergence  will  be  obtained  when  the  optimal  racial  meshes  with  optimal 
number  of  elements  are  used  and  p  ->  oo? 

Theorem  2.4.4  [36].  There  exists  sequence  of  meshes  A  (dependent  on 
p)  such  that  for  uniform  degrees  of  elements  we  have 

1  -4V(a-l/2)N  /e 

This  error  is  obtainable  by  using  the  proper  sequence  of  radical  meshes  with 

p  w  v'4(a-l/2)N/e  as  N  oo,  0  =  .  Let  us  summarize  the  theorem  in 

a-1/2 

Table  2.4.1.  For  r  =  ( l-v'q)/(  1+Vq)  we  have 

Z 


(2.4.15) 


1  K^(a-1/2)N 


/iT 


(2.4. 14) 


■Z“E 


TABLE  4.1.  The  error  for  the  h-p  version. 


METHOD 

cr 

s 

K 

0- 

Geometric 

mesh 

nonuniform 
p-dlstribution 
Theorem  2.4.2 

q 

tn  r 

VZ  £n  q  £n  r 

0 

1/2 

0.3932 (a-1/2) 

1 . 5632 

0 

(v^-1)  ^ 

2a  -  1 

1 . 7627 

0 

Geometric 
mesh 
uniform 
p-distribution 
Theorem  2.4.3 

q 

Vtn  q  tn  r 

min(a, 2a-l ) 

1/2 

0.3932 (a-1/2) 

1.1054 

min(a, 2a-l ) 

2a  -  1 

1 . 2464 

min(a, 2a-l ) 

Radical  mesh 
q  and  s 
asymptotic 
Theorem  2.4.4 

1.4715 

a  -  1/2 
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Let  us  now  briefly  analyze  Table  4.1.  The  h-p  version  leads  to  the 

...  .  ^  -KVTo-T/ZTN  r-  i.1,  T  ^ 

exponential  rate  of  convergence  e  For  the  uniform  p  and 

optimal  geometrical  mesh  k  =  1.24  while  for  the  optimal  nonuniform  mesh  we 
get  K  =  1 . 76  =  1 . 24  .  Analogous  results  hold  in  2  dimensions  when  the 

uniform  degrees  leads  to  the  exponents  which  is  by  factor  smaller  than 

the  exponent  for  the  nonuniform  degree  distribution  on  a  geometric  mesh.  For 
uniform  degrees,  a  sequence  of  the  meshes  with  the  exponent  1.47  can  be 
constructed.  The  shown  results  indicate  roughly  what  kind  of  meshes  should  be 
considered  for  the  mesh  design  for  the  p-version  . 

The  mentioned  results  have  to  be  interpreted  in  the  light  of  adaptive 
approaches  when  applied  for  practical  computations. 

Once  more  two  strategies  of  the  mesh  generation  could  be  used.  The 
local  and  global  one  as  discussed  in  Section  2.2.  In  the  local  approach  the 
main  difficulty  is  to  decide  whether  change  of  the  degree  or  refinement  (or 
derefinement)  of  an  element  has  to  be  made.  See,  e.g.  [36]  [48]  [50]  for  a 
discussion.  In  the  global  strategy,  the  main  difficulty  is  to  extract  from  the 
solution  the  information  allowing  to  design  global  mesh  and  global  degree 
distribution.  The  adaptive  approach  based  on  p  uniform  which  is  made  as  a 
sequence  of  meshes  for  increasing  uniform  degrees  (analogously  as  in  Theorem 
2.4.4)  seems  too  impractical. 

The  h-p  version  was  theoretically  analyzed  in  the  two  dimensional  cases 
in  [4] [6] [8] [9] [16] [11] [15] [21] [36] [37] [38] [48] [55] [59] .  See  also  the  survey 
[20]. 

3.  THE  FINITE  ELEMENT  METHOD  IN  TWO  DIMENSIONS 
3.1.  Introduction 

In  this  section  we  will  discuss  the  performance  of  the  p  and  h-p  version 
for  the  elliptic  boundary  value  problems  with  piecewise  analytic  data.  We 
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will  consider  polygonal  domains  with  straight  or  curved  sides,  the  piecewise 
analytic  coefficients  of  the  operator,  piecewise  analytic  right  hand  side  and 
boundary  conditions  which  type  could  change.  This  class  of  problems  is 
typical  in  engineering 

It  is  well  known  that  the  solution  has  singular  behavior  in  the  corners 

of  the  domain,  places  where  the  boundary  condition  abruptly  changes,  etc. 

See  [27] [37] [42] .  Although  the  solution  behaves  in  this  areas  similarly  as 

the  function  (2.1.2)  in  Section  2.1,  the  situation  here  is  much  complex. 

Mathematical  tools  of  accurate  description  of  the  regularity  of  the  solution 

k  t 

are  needed.  Such  tool  is  the  weighted  Sobolev  space  H-’  (fl)  and  the 

P 

I 

countably  normed  space  The  use  of  these  spaces  allows  to  prove  very 

P 

accurate  estimates  of  the  error  of  the  p  and  h-p  version  of  the  finite 
element  method  which  are  very  analogous  to  the  one  cited  in  Section  2, 
although  not  so  detailed.  To  asses  the  detailed  applicability  of  the 
theoretical  results  for  practical  computation,  we  will  present  various 
numerical  examples  and  tests. 

The  h-p  and  p  versions  in  practical  environment  have  to  be  understood  in 
a  context  of  adaptive  procedure  with  an  error  estimator.  We  will  address 
these  aspects  too. 

We  will  not  elaborate  here  on  the  h-version  because  this  is  a  classical 
approach  and  many  results  are  available. 


k  £  I 

3.2.  The  spaces  (n)  and  and  the  model  problem 


Let  n  c  Ik  be  a  polygon  with  vertices  A  and  (open)  edges  F,  = 

M  ^ 

^  s  i  2!  M,  ^^+1  “  9£2  =  U  Fj^ .  By  we  denote  the 

interior  angle  of  A^  and  assume  0  <  s  2n.  Let  r^  (x)  =  dist(x,A^)  and 

B  =  O,  ,...,|3u)  be  M-tuple  of  real  numbers  e  (0,1).  Further,  for  any 
^  ”  M  k^B^  ^ 

integer  k  we  define  ♦^^^^(x)  =  ]  \  r^  (x). 
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Let  H  (Q)  be  the  standard  Sobolev  space.  Then  for  k  ^  ^  0, 

k  ^ 

integers  we  define  the  space  H^’  (Q)  as  the  completion  of  the  set  of  all 
infinitely  differentiable  functions  under  the  norm 


I  Ir  i 

H^’  (fl) 
P 


1^11  £_1  ^ 


Jsla  sk 


OS  a  sk 


We  further  introduce  the  countably  normed  space 

B^(£2)  =  {u|u  €  Vks£.  I V  |  a| -£  (Cl)  '  '  • 

lal  =  k,k+l .  constants  Cal  and  dal  independent  of  k> 


Remark  3.2.1.  The  intuition  behind  the  definition  of  the  spaces 
V  t  t 

Hr’  (n)  and  Bo(n)  is  to  consider  function  u  =  r  in  the  neighborhood  of 

p  p 

the  corners  and  achieve  that  it  naturally  belongs  to  these  spaces.  The 
solution  u  behaves  in  the  neighborhood  of  the  corners  similarly  as  this 
function.  See  [7] [ 1 1 ] [ 12] [37] . 

Let  us  consider  the  following  model  problem. 


(3.2.1) 


-Au  =  f  on  n 


“'Fo  '  ■  Slfp 

"lr„  '  %  =  Slr^ 


where 


U  r  r 

IsE 


U  r,  ;  V  *  0  is  a  subset  of  M  =  (1,2, ... ,M} 
i€A'  ^ 


and  A  =  M  -  D.  Let  =  {u|u  €  H^(Q),u  =  0  on  r^}.  The  weak  form  of 

(3.2.1)  is  u  e  H^,  u  =  Gj^  on  and  such  that  for  any  v  €  Hp(n) 
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(3.2.2) 
with 

(3.2.3) 


B(u, v) 


r  fv  dx  + 
•^Q 


G„v  ds 


F(v) 


B(u,v)  =  r  7u7v  dx 


Similarly,  as  in  one  dimension,  we  define  the  energy  norm 


(3.24) 


||u||^  =  (B(u,u))^^^ 


which  is  equivalent  to  1*11  on  H^(Q). 

H^Q) 

Consider  now  the  domain  fl  shown  in  Fig.  3.2.1.  Then  we  have 
Theorem  3.2.  t  171.  If  feH'‘’°(a).  G.  e  G..  e  ■  '  (B) 

I  “  @  "  g 

with  €  (0,1),  1  s  i  s  M,  then  problem  (3.2.1)  has  unique  solution  u  e 

(^■^^’^(n)  with  =  0,  if  >  1  -  K.  or  3.  €  (l-ic,,l)  if  p.  ^  1  - 

K,  ,  where  k,  =  —  if  r,,r.  ,  c  Fr,  (or  F.,)  and  k,  =  ^  otherwise, 
i  i  w.  i  1+1  D  N  i  2u>, 


Fig.  3.2.1.  Polygonal  domain  £2. 

Further,  if  f  e  B°(n).  G„  e  B^(Q)  and  G„  e  s'ln),  then  u  e  B^ln). 

B  “  B  "  S  ® 

Remark  3.2.2.  The  weight  depends  on  and  which/is  related 

to  the  geometry  of  Q.  Hence,  also  if  f,  and  Gjj  are  analytic,  the 
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solution  still  may  have  singular  behavior  because  of  the  unsmootheness  of  the 
boundary  f  of  □ 

Remark  3. 2. 3.  If  £2  is  a  polygon  with  the  analytic  curvilinear  edges, 

the  solution  u  of  (3.2.1)  will  be  in  (£2),  c  >  0  (instead  of  B^(£2)) 

B+e  B 

with  c  >  0  arbitrary  (see  [7]) 


Remark  3.2.4.  Theorem  3.2.1  although  formulated  only  for  the  problem 
(3  2.1)  holds  also  for  other  problems  as  for  the  elasticity  problems  where  the 
values  of  B^^  and  Bj^  are  properly  adjusted  (i.e.,  they  are  not  the  same  as 
in  Theorem  3.2.1).  Analogous  results  are  also  valid  for  the  eigenvalue 
problems,  etc.  For  more,  see  (121. 

In  the  next  section  we  will  address  the  problem  (3.2.1)  and  analogous 
problems  for  the  elasticity  and  will  always  assume  that  Gj^  =  0. 


3.3.  The  h-p  version 

The  h-p  version  simultaneously  refines  the  mesh  and  increases  the 

degrees  of  the  elements.  In  this  section  we  assume  that  the  solution  u  of 

the  problem  has  singularity  only  at  one  vertex  of  the  domain.  Let  us  define 

in  the  neighborhood  w  of  the  vertex  the  geometrical  mesh  £2^  =  1  £  i 

sm,  lsjij(i)>  where  £2^^^  are  the  elements,  the  triangles  or 

quadrilaterals  (straight  or  curved).  Example  of  a  sequence  of  geometrical 

meshes  is  shown  in  Fig.  3.3.2.  The  elements  £2.  ,  satisfy  the  usual 

ij 

conditions  preventing  their  degeneration  and  other  standard  conditions. 

We  denoted  by  a  the  mesh  factor  (denoted  in  Section  2  by  q)  and  by 


m  the  number  of  layers.  Elements  £2^^^,  J  =  1,2 . J(i)  are  located  in  the 

(i-l)-th  layer.  Denote  by  d.  .  the  distance  between  the  element  £2.  .  and 

ij  ij 

the  vertex  and  by  h.  .  and  h.  .  the  maximal  and  minimal  side  of  the  element 

ij  “iJ 


£2.  We  will  assume  that  d..  «h.,  ~h..  =o- 

iJ-  ij  ij  -ij 


m-i 


l<ism,  lij2j(i). 
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djj  =  0,  ~  IsjsJ(l).  Further,  we  will  call  the  vector 

p  =  ^Pij’  l-i-m,  l£j£  J(i)},  Pj^j  >  1  integer  the  degree  vector.  We 

defined  only  the  elements  in  the  m  layers  near  the  vertex  covering  w. 

On  n-w  we  complement  the  mesh  {n"*!  by  the  elements  {n  >  =  Let  IM.  .} 

<r  r  O'  ij 

be  the  mappings  of  the  standard  elements  S  onto  resp.  onto  We 

will  assume  that  these  mappings  satisfy  the  usual  conditions  of  the  finite 
element  method.  The  finite  element  space  S^’^(n)  is  then  defined  by 


5^(0^)  =  =  *(M  "(x.y)), 


-1 


is  a  polynomial  of  degree  p^^j  on  S> 

and  analogously  on  {0^}  with  the  degree  p^  =  max  S^’^(n)  =  © 

sP,c^mjj  ^  and  =  S^’^CQ)  n  H^(n).  As  usual,  the  finite 

element  solution  u^  e  §^’^(0)  satisfy 

B(u  ,v)  =  B  (u.v)  ^  F^v),  Vv  € 

s 

and  we  set  Theorem  3.3.1  [8].  Let  Q*”  be  the  geometrical  mesh  with  the  mesh 

-  {T 

factor  <r  €  (0,1)  and  let  p  >  0  be  a  degree  factor  such  that  p_  =  p^  =  1 
+  [p(i-l)]  (linear  distribution  )  or  p  =  1  +  [pmj  (uniform  degree 

distribution),  and  p  =  max  p,  .  on  .  If  u  e  B^(J2)  n  Hi(n),  is  the 

1.1  O’  p  u 


solution  of  the  problem,  then 
(3.3.1) 


||u-Ug||^  - 


-bN 


1/3 


where  N  is  the  number  of  degrees  of  freedom  (i.e.  dimension  of  S^’^(n)) 
and  constant  C  and  b  are  independent  of  N.  We  see  an  analogous  result  as 
in  Theorem  2.4.2  but  without  exact  specification  of  the  constants  C  and  b. 

Let  us  illustrate  Theorem  3.3.1  by  a  numerical  example.  Consider  the 
elasticity  problem  on  a  cracked  domain  shown  in  Fig.  3.3.1.  We  will  assume 
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that  the  material  is  homogeneous  and  Poisson  ratio  v  =  0.3.  The  exact 
solution  is  the  first  stress  intensity  mode  and  tractions  are  prescribed  on 
r  6  Q.  Because  of  the  symmetry,  we  consider  only  the  upper  half  of  the 
domain.  The  sequence  of  the  used  meshes  with  o-  =  0.15  and  n  =  m  +  1 

is  shown  in  Fig.  3.3.2.  The  uniform  degree  distribution  was  used  with  p  = 

1  +  [pm],  p  =  1.  The  computation  has  been  made  by  the  program  MSC/PROBE.  In 
Table  3.3.1  we  show  the  basic  results  and  the  relative  error  ||e 

100%  with  b  and  C  in  (3.3.1).  By  1|Uq11£  we  denoted  the  energy  norm  of 


the  exact  solution.  In  Fig.  3.3.3  we  show  the  error  as  function  of  N.  We 
present  here  also  the  error  for  the  p-version.  The  error  of  the  h-p  version 
(5=1)  is  shown  by  the  solid  line.  We  see  that  the  line  is  straight  as 
expected  from  the  theory.  The  values  of  constants  b  and  C  cannot  be  at 


Fig.  3.3.1.  Cracked  panel. 


Figure  3.3.2.  Geometric  mesh  A^.  isn^b. 
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present  theoretically  determined.  In  Fig.  3.3.4  we  show  the  error  of  the  h-p 
version  for  c  =  1  and  various  <r.  We  see  that  the  value  cr  =  0.15  is  nearly 
optimal . 

TABLE  3.3.1 


Performance  of  the  h-p  version  on  mesh  A^,  1  s  n  s  6,  o'  =  0.15,  p  =  1 . 


Mesh 

P 

N 

^1/3 

■n 

b 

1 

9 

60.92 

0.741 

^2 

2 

48 

3.63 

20.23 

0.740 

2.303 

^3 

3 

121 

4.95 

7.61 

0.776 

4 

256 

6.35 

2.57 

0.720 

1.810 

^5 

5 

477 

7.82 

0.670 

1.683 

6 

808 

9.31 

0.33 

0.670 

1.688 

23456789  10  n 


Fig.  3.3.3.  Performance  of  the  p  and  h-p  version 
on  Mesh  A,lsns6,  o-  =  0.15. 
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Fig.  3.3.4.  Performance  of  geometric  Mesh 
with  <r  =  0.15,  0.08  and  0.5. 


We  note  that  the  results  are  analogous  to  those  in  Section  2.  If  the 
degrees  of  the  element  are  linearly  distributed  we  would  get  b'  ~  v  3  b  in 
(3.3.1)  and  Table  3.3.1. 

The  adaptive  procedure  consists  of  the  design  of  the  estimator  E  and 
the  adaptive  principles.  The  error  estimator  can  be  based  on  the 
extrapolation  technique.  In  Table  3.3.2  we  show  the  exact  error  ||e||^  ,  the 
computed  error  estimator  ©  =  ||e]|£  and  the  relative  error  of  the  estimator  in 
%.  We  see  very  high  reliability  of  the  estimator.  The  adaptive  approach  for 
the  h-p  version  is  a  complex  one.  Here  once  more  the  local  and  glv->bal 
adaptive  approach,  as  explained  in  Section  2,  can  be  used.  In  contrast  to  the 
one  dimensional  case  transitional  elements  between  the  element  of  different 
degrees  are  needed.  In  addition,  various  shape  functions  can  be  adaptively 
chosen  so  that  elements  are  not  of  a  full  degree.  Further,  different  type  of 
elements,  triangular  or  quadrilateral  can  be  used  and  the  degree  of  element 
can  also  be  understood  as  a  total  or  tensor  product  type. 
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TABLE  3.3.2 


Error  estimator  of  the  h-p  version  on  Mesh  A^,  1  s  n  £  6,  <r  =  0.15,  p  =  1. 


mm 

I'Ie.r''- 

(l|ellE-|ellE)/||e||^% 

2.9596E-1 

2.9662E-1 

60.83 

60.92 

0.2189 

9.8774E-2 

9.8511E-2 

20.28 

20.23 

-0.2669 

3,7033E-2 

3.7055E-2 

7.606 

7.611 

0.0606 

1 . 2489E-2 

1 . 2500E-2 

2.565 

2.567 

0.0926 

4.3359E-3 

4.3691E-3 

0.891 

0.897 

0.6689 

The  implementation  influences  very  heavily  the  adaptive  process.  For  some 
results,  we  refer  here  to  [2] [3] .  It  seems  once  more  that  the  h-p  adaptive 
approach  based  on  the  adaptively  constructed  meshes  for  the  h-version  and 
various  p  (analog  to  the  case  addressed  in  Theorem  2.4.4)  is  practically  not 
usable.  Once  more  we  see  that  the  theoretical  results  could  give  a  very  good 
insight  in  the  problem  of  the  design  of  an  adaptive  procedure.  The  first 
theoretical  analysis  of  the  h-p  version  in  conjunction  with  the  regularity  of 
the  solution  described  by  countably  normed  spaces  appeared  in  [38],  For 
more  results,  we  refer  to  [8] [9] [10] [11] [12]  and  the  survey  [20]  where 
additional  references  are  available. 

3.4.  The  p  version 

In  the  p-version  the  mesh  is  fixed  and  the  degrees  of  the  elements  are 

increased  uniformly  or  nonuniformely.  If  the  domain  Q  has  corners  and  the 

solution  is  singular  in  the  areas  of  these  corners,  the  rate  of  convergence  of 

the  method  is  algebraic  (see  [13]  [20]  [21]  [22]).  This  is  quite  analogous 

2  -2(1-6) 

behavior  results  as  in  one  dimension.  If  u  e  ^^(n)  then  ||e||£  -  CN  . 

The  rate  is  twice  as  large  as  in  the  case  of  the  h-version  with  a  quasiuniform 
mesh.  In  the  case  that  the  solution  is  analytic  on  Q,  we  have  ||e||^  s 
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-bN 


1/2 


-bp 


e  «  Ce  .  As  in  one  dimension,  the  accuracy  is  for  large  p  governed 

by  the  error  on  the  element  containing  the  vertex  of  the  domain,  but  for 
smaller  p  and  the  refined  mesh,  the  error  could  be  governed  by  the  elements 
far  from  the  singularity. 

As  was  said  in  Section  2,  the  user  designs  the  mesh  which  is  refined  in 

the  neighborhood  of  the  corners  and  other  critical  places  and  the  degree  of 

the  elements  is  determined  adaptively.  The  mesh  which  typically  has  not  too 

many  elements  is  designed  by  the  experience  and  simple  rules.  Another  way  is 

to  use  a  crude  mesh  constructed  directly  by  the  user  or  adaptively  by  the 

h-method  with  small  p,  say  3.  Then  this  mesh  could  be  further  refined  by 

the  users  experience  (see  e.g.  [59])  or  by  an  expert  system  (see,  e.g.  [15]). 

We  underline  that  the  degrees  of  the  elements  are  adaptively  determined  so 

that  reliable  results  are  obtained  for  any  choice  of  the  meshes;  of  course,  a 

"bad"  mesh  will  lead  to  a  higher  computer  cost. 

We  will  assume  that  the  constructed  meshes  are  of  geometric  type  in  the 

neighborhood  of  the  vertices  with  m^  layers,  i.e.,  we  will  consider  the 
">0 

mesh  with  m^  fixed.  From  the  analysis  of  [8]  we  have 


(3.4.1) 


Me  ^  c 


r 

0" 


where  0  <  Q  <  1  and  y  <  0. 

The  first  term  characterizes  the  error  in  the  elements  closest  to  the 
singularities  and  the  second  term  all  others.  For  large  p  (as  in  one 
dimension),  the  first  term  dominates.  For  small  p  the  second  one  dominates 
and  the  rate  is  exponential.  Hence,  the  p-version  has  two  phases,  asymptotic 
and  pre-asymptotic.  The  preasymptotic  phase  is  until  the  elements  closest  to 
the  singularity  is  not  governing  the  accuracy.  Hence,  the  user  tends  to 
design  such  a  mesh  for  which  the  desired  accuracy  in  both  terms  in  (3.4.1)  is 
balanced. 

Let  us  present  the  computational  results  for  the  same  elasticity  problem 
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as  in  Section  3.3.3.  In  Fig.  3.3.3  we  already  have  shown  the  results  of  the 
p-version.  In  Fig.  3.4.1  we  show  the  relative  error  in  the  In  x  in  scale. 

N 


Fig.  3.4.1.  The  error  of  the  p-version. 

We  see  that  for  a  fixed  mesh  the  accuracy  curve  has  a  S  shape  with  the 
pre-asymptotic  and  asymptotic  phase  and  inflection  pint  roughly  when  the  error 
of  the  elements  closest  to  the  vertex  become  to  be  dominant. 

Let  us  now  consider  the  problem  for  the  Laplace  equation  on  the  domain 

1/2 

shown  in  Fig.  3.3.1  and  the  exact  solution  u  =  r  cos  0/2,  where  (r,0) 

are  the  polar  coordinates.  Because  of  the  symmetry  only  one  half  of  the 

domain  will  be  considred.  The  boundary  condition  is  of  Neumann  type.  On  the 
side  {-1  <  X  <  0,  y  =  0}  zero  Dirichlet  condition  is  prescribed.  The  used 

mesh  Is  shown  in  Fig.  3.3.2  except  the  two  square  elements  containing 

singularity  are  replaced  by  four  equal  triangular  elements.  We  use  <r  =  0.13 
and  <r  =  0.50. 

In  Figs.  3.4.2  and  3.4.3  we  show  the  relative  error  ||6i£  p  ’/» 
function  of  N  in  the  in  x  in  scale  for  different  p  and  various  number  of 
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error  is  not  decreasing  for  the  increased  number  of  layers.  This,  as  in  one 
dimension,  depends  on  the  degree  p  of  the  elements.  We  see  that  the  number 
of  3-4  layers  and  o*  =  0.13  is  optimal  and  gives  very  good  performance.  For 
p  increasing  the  case  cr  =  0.5  is  also  leading  to  the  convergent  results  but 
less  efficiently.  So  we  see  that  overrefinement  is  desirable. 

In  Tables  3.4.3,  3.4.4  and  3.4.5  we  compare  the  performance  of  the 
meshes  for  different  ranges  of  accuracy. 

Tables  3.4.3  -  3.4.5  show  well  the  influence  of  the  number  of  layers  on 
the  accuracy.  We  emphasize  that  the  element  degree  selection  leading  to  the 
desired  accuracy  is  made  adaptively.  In  practice  m  »  3  is  usually  used  and 
here  o-  «  0.15  is  a  good  choice  for  the  accuracies  of  1-3%.  To  show  it,  we 
present  in  Table  3.4.6  the  achieved  accuracy  for  m  =  3  as  function  of  p 
and  N  for  <r  =  0.13  and  <r  =  0.50. 

The  computations  results  we  have  presented  are  in  very  good  agreement 
with  the  theoretical  analysis. 
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TABLE  3.4.1.  Performance  of  the  p  version  on  geometric  mesh  with  (t  =  0. 13 


Mesh 

Number  of 

mgm 

p  =  2 

A 

n 

layer,  m 

N 

N 

N 

m 

N 

I'Ie,/- 

A 

2 

1 

8 

29.  13 

24 

13.03 

44 

8.33 

72 

6.21 

A 

3 

2 

12 

23.53 

36 

8.20 

64 

4.00 

104 

2.49 

H 

3 

16 

22.48 

48 

7.28 

84 

2.98 

A 

5 

4 

20 

22.29 

60 

7.15 

104 

2.82 

168 

A 

6, 

5 

24 

22.72 

72 

7. 13 

124 

2.80 

200 

BB 

B 

6 

28 

22.26 

84 

7. 13 

144 

2.80 

232 

BQI 

A 

8 

7 

32 

22.26 

96 

7. 13 

164 

2.80 

264 

BB 

A 

9 

8 

36 

22.26 

108 

7. 13 

184 

2.80 

296 

B 

Mesh 

Number  of 

p  =  5 

■BbB 

p  =  7 

p  =  8 

A 

n 

layer,  m 

N 

m 

N 

m 

N 

A 

2 

1 

108 

5.03 

152 

4.24 

204 

3.67 

264 

Bl 

A 

3 

2 

156 

1.87 

220 

1.54 

296 

1.33 

384 

1.17 

B 

3 

204 

0.81 

288 

0.59 

388 

0.48 

504 

0.42 

A 

5 

4 

252 

0.54 

480 

0. 19 

624 

0.  15 

A 

6 

5 

300 

0.547 

572 

0. 11 

744 

0.068 

B 

6 

348 

0.49 

492 

0.21 

664 

0.094 

864 

0.  046 

A 

8 

7 

396 

0.49 

560 

0.21 

756 

0.093 

984 

0.041 

A 

9 

8 

444 

0.49 

268 

0.21 

847 

0.092 

1104 

0.040 
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TABLE  3.4.2.  Performance  of  the  p  version  on  geometric  mesh  with  <r  =  0.50 


Mesh 

A 

n 

Number  of 
layer,  m 

p  =  1 

p  =  2 

IHBiH 

p  =  4 

N 

N 

N 

BB 

BB 

A 

2 

1 

8 

36.01 

24 

21.40 

44 

15.34 

72 

11.96 

A 

3 

2 

12 

27.45 

36 

15.30 

64 

10.89 

104 

8.48 

n 

3 

16 

21.60 

48 

10.95 

84 

7.72 

136 

6.00 

A 

5 

4 

20 

17.85 

60 

7.88 

104 

5.47 

168 

mm 

A 

6 

5 

24 

15.58 

72 

7.76 

124 

3.88 

200 

3.00 

B 

6 

28 

14.30 

84 

4.31 

144 

2.75 

232 

2.  12 

A 

8 

7 

32 

13.60 

96 

3.36 

164 

1.96 

264 

1.50 

A 

9 

8 

36 

13.23 

108 

2.76 

184 

1.41 

296 

IBI 

A 

16 

9 

40 

12.95 

120 

2.41 

204 

1.03 

328 

B 

B 

10 

132 

2.21 

224 

0.77 

360 

0.53 

A 

12 

12.90 

144 

2.11 

244 

0.60 

392 

0.38 

1? 

56 

12.86 

168 

=== 

2.02 

284 

0.44 

456 

0.20 
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TABLE  3.4.2.  Performance  of  the  p  version  on  geometric  mesh  with  <r  =  0.50. 

(Continued) 


Number  of 
layer,  m 


m 


p  =  5 


N  N 


p  =  6 


ellE.n/-  N 


p  =  8 


108  9.80  152  8.30  204  7.20  264  6.35 


156  6.94  220  5.87  296  5.09  384 


204  4.91  288  4.16  388  3.60  504  3.18 


252  3.47  356  2.29  480  2.55  624  2.25 


300  2.45  424  2.08  572  1.80  744  1.59 


348  1.74  492  1.47  664  1.27  864  1.126 


396  1,23  560  1.04  756  0.90  984 


444  0.87  628  0.73  848  0.64  1104  0.56 


492  0.61  696  0.52  940  0.45  1224  0.40 


540  0.43  764  0.37  1032  0.32  1344  0.28 


588  0.31  832  0.26  1124  0.23  1464 


588  0.15  968  0.13  1308  0.11  1704 


TABLE  3.4.3.  Comparison  of  performance  of  the  p-version  on  n 


for  the  accuracy  1%. 
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TABLE  3.4.4.  Comparison  of  the  p-version  performance  on  n"*  for  the  accuracy  3% 


p 

O'  =  0.  13 

O'  =  0.5 

m 

N 

m 

N 

2 

8 

108 

2.76 

3 

3 

84 

2.98 

6 

144 

2.75 

4 

2 

104 

2825 

5 

200 

3.00 

5 

2 

156 

1.87 

5 

300 

2.45 

6 

2 

220 

1.54 

4 

356 

2.94 

7 

2 

296 

1.33 

4 

480 

2.55 

8 

1 

264 

3.24 

3 

504 

3.  18 

TABLE  3.4.5.  Comparison  of  the  p-version  performance  on  Q*"  for  the 
accuracy  12%  . 


P 

O'  =  0.13 

O'  =  0.  15 

m 

N 

I'Ie.r-''- 

m 

N 

||e|l£_R% 

1 

10 

44 

12.95 

2 

1 

24 

13.03 

3 

48 

10.95 

3 

1 

44 

8.33 

2 

64 

11.89 

4 

1 

72 

6.21 

1 

72 

11.96 

5 

1 

108 

5.03 

1 

108 

9.80 

6 

1 

152 

4.24 

1 

152 

8.30 

i|Bil 

1 

204 

3.67 

1 

204 

7.20 

1  ^ 

1 

264 

3.24 

1 

264 

6.35 

TABLE  3.4.6.  Comparison  of  performance  of  the 
p  version  on  n  . 

O' 


p 

N 

nn 

1 

1 

■OB 

1 

30 

22.47 

21.10 

2 

60 

7.28 

10.95 

3 

104 

2.97 

7.72 

4 

168 

1.37 

6.00 

5 

252 

0.76 

wmm 

6 

356 

0.51 

4.  16 

7 

480 

0.38 

3.60 

8 

674 

0.30 

3.  18 

TABLE  3.4.7.  Performance  of  Error  Estimator  on  Mesh  n"' 
<r  =  0.  13,  m  =  2. 


P 

N 

■BBI 

0 

1 

12 

23.52 

23.56 

1.00 

2 

36 

8. 19 

8.26 

0.999 

3 

64 

3.97 

4.00 

0.992 

4 

104 

2.45 

2.49 

0.984 

5 

156 

1.82 

1.87 

0.973 

6 

220 

1.47 

1.34 

0.954 

7 

296 

1.25 

1.32 

0.947 

8 

384 

1.08 

1.  17 

0.923 

As  we  said,  the  essential  part  of  the  code  is  an  error  estimator.  For 
the  p-version  the  estimator  based  on  the  extrapolation  is  very  effective.  The 
quality  of  the  estimator  can  be  judged  by  its  effectivity  index  0 

g  _  estimate  _  S 
true  error  ~  |e  || 

which  can  be  computed  by  in  benchmark  computations.  We  show  in  Table  3.4.7 
the  quality  of  the  estmator  for  on  example  when  the  mesh  (i.e., 

m  =  2)  is  used  for  a-  =  0.13. 

3.5.  The  adaptive  p  and  h-p  versions 

The  aim  of  the  adaptive  approach  is  to  obtain  (in  a  most  effective  way) 
the  solution  in  the  a-priori  given  range  of  accuracy  with  a  reliable  error 
estimate.  The  effectiveness  is  meant  in  general  way,  i.e.,  includes  the 
computer  and  manpower  cost.  The  adaptive  approach  is  very  much  influenced  by 
the  implementation  aspects.  The  code  based  on  the  h,  p  and  h-p  version  have 
very  different  characters.  Here  we  will  concentrate  on  the  p  version  and  the 
error  measure  of  the  energy  norm. 

1)  The  mesh  generation  is  for  complex  geometry  very  laborious.  Various 
mesh  generators  commercial  and  not  commercial  are  available  today.  They  focus 
primarily  on  the  handling  of  the  geometry.  In  the  case  of  the  h-version  (in 
two  dimensions)  various  adaptive  codes  are  available.  Nevertheless,  these 
meshes  are  not  appropriate  for  the  p-version  (although,  of  course,  they  can  be 
used  too). 

2)  The  mesh  for  the  p-version  is  characterized  by  large  elements  with  a 
(geometric)  refinement  in  the  neighborhood  of  singularities. 

3)  The  mesh  in  the  neighborhood  of  the  singularity  has  to  be  rather 
overrefined  than  underrefined.  The  underrefinement  can  be  more  "risky".  The 
reasonable  overref inement  is  not  costly.  The  underref inement  can  lead  to 
pollution  effects.  For  the  stronger  singularity  more  layers  have  to  be  used. 
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Usually  we  will  recommend  2-4  layers..  The  factor  <r  of  the  geometrical  mesh 
in  the  neighborhood  of  the  corners  can  be  different,  nevertheless,  we 
recommend,  for  the  reasons  explained  in  previous  sections,  to  use  c  ~  0.15. 

4)  The  mesh  generator  has  to  have  special  features  to  lead  to  the 
meshes  which  are  advantageous  for  the  p-version.  They  can  be  based  on  various 
principles.  For  example,  on  the  design  of  the  crude  mesh  and  special  "layer" 
refinement  in  the  neighborhood  of  the  critical  places  possibly  by  help  of  an 
expert  system  (see,  e.g.  [15]),  etc. 

5)  If  the  information  about  the  solution  in  some  point  the  neighborhood 
of  the  corner  is  needed,  then  these  points  have  to  be  separated  from  the 
corner  by  2  layers.  The  stresses  in  the  area  very  close  to  the  corners  have 
to  be  computed  by  stress  intensity  factors. 

6.  The  p-version  could  suffer  by  an  oscillation  of  the  solution 
if  the  mesh  is  improperly  designed.  This  could  happen  in  the  elements  which 
contain  the  singularity.  This  element  has  to  be  small  and  direct  stress 
computation  in  this  elements  should  not  be  used.  We  note  that  the  energy 
error  converges  also  in  this  element.  Nevertheless,  experience  backed  by  the 
theory  shows  that  for  proper  mesh  no  oscillation  occurs. 

7)  The  error  indicators  and  estimators  based  on  the  comparison  of  data 
for  various  consecutive  p  performs  very  well,  not  only  for  the  energy  norm 
accuracy  measure.  This  is  essential  because  the  energy  norm  measure  is 
usually  not  the  most  important  one  and  possibly  could  be  misleading  when  other 
accuracy  measure  is  of  interest. 

8)  The  error  indicator  can  be  used  elementwise.  This  is  essential  if 
the  accuracy  will  not  be  achieved  in  the  range  of  admissible  range  of  degrees 
p.  Then  a  new  mesh  is  created  by  refinement  of  the  elements  with  largest 
error  indicators. 

9)  So  far  we  mentioned  only  the  p  adaptivity  with  uniform  p.  The 


38 


adaptive  approach  can  be  based  on  adaptive  additions  of  various  shape 
functions  so  that  the  notion  of  degrees  is  loosing  its  meaning. 

10)  The  experience  holds  also  for  the  p-version  that  quadrilateral 
elements  are  better  than  the  triangular  ones. 

11)  One  of  the  main  advantage  of  the  p-version  is  relatively  easy 
implementation  and  possibility  of  the  reliable  error  control  for  any  data  of 
interest  by  comparing  the  computed  data  for  various  p. 

12)  The  implementation  and  adaptive  strategy  for  the  h-p  version  is  much 
more  complicated  but  could  be  made,  see,  e.g.  PHLEX  program  [49]  [50].  The 
possible  plausible  approach  to  construct  the  h-p  version  by  optimal  meshes  for 
the  h-version  for  increasing  p  cannot  be  recommended  for  practical  reasons. 

13)  For  additional  aspects  of  the  p-version  for  the  analysis  of 
complicated  engineering  problems  by  the  STRIPE  program,  we  refer  to  [2]. 

4.  CONCLUSIONS 

We  presented  basic  ideas  behind  the  h,  p  and  h-p  versions  and 
characteristic  theoretical  results.  The  comparison  of  these  approaches  is  a 
complex  task  dependent  on  the  criteria  used.  It  reflects  the  implementation 
aspects,  computer  hardware,  etc.  By  our  opinion  one  of  the  basic  comparison 
criterion  has  to  be  the  possibility  of  accurate  and  reliable  quantitative 
assessment  of  the  accuracy  of  any  data  of  engineering  interest.  It  seems  that 
the  adaptive  p-version  (on  a  properly  designed  mesh)  is  one  of  the  best 
approach  for  the  practical  engineering  computations  of  problems  in  structural 
mechanics.  (See  also  [2].)  Of  course,  changing  criteria  of  comparison  or  a 
class  of  the  problem  could  lead  to  different  conclusions. 


39 


REFERENCES 


[1]  Ainsworth,  M. ,  Oden,  J.  T.  [1991];  A  Unified  Approach  to  A  Posteriori 
Error  Estimation  Based  on  Element  Residual  Methods.  Submitted. 

[2]  Andersson,  B. ,  Babuska,  I.,  Petersdorff,  v.T.  [1992]:  Reliable  Stress 
and  Fracture  Mechanics  Analysis  of  Complex  Aircraft  Component  Using  a 
h-p  Version  of  FEM.  To  appear. 

[3]  Andersson,  B. ,  Babuska,  I..  Petersdorff.  v.T.  [1992]:  Computation  of 
the  Vertex  Singularity  Factors  for  Laplace  Equation  in  3  Dimensions. 

[4]  Babuska,  I.,  Door,  M.  [1981]:  Error  Estimate  for  the  Combined  h  and  p 
Versions  of  the  Finite  Element  Method,  Numer.  Math.  37,  257-277. 

[5]  Babuska,  I.,  Duran,  R. ,  Rodriguez,  R.  [1992]:  Analysis  of  the 
Efficiency  of  an  a-posteriori  Error  Estimator  for  Linear  Triangular 
Finite  Elements,  SIAM  J.  Numer.  Anal.  To  appear. 

[6]  Babuska,  I..  Gui,  W.  [1986]:  Basic  Principles  of  Feedback  and  Adaptive 
Approaches  in  the  Finite  Element  Method.  Comput .  Methods  Appl.  Mech. 
Engrg.  55,  27-42. 

[7]  Babuska,  I.,  Guo,  B.Q.  [1988]:  Regularity  of  the  Solutions  of  Elliptic 
Problem  with  Piecewise  Analytic  data.  Part  I:  Boundary  Value  Problems 
for  Linear  Elliptic  Equation  of  Second  Order,  SIAM  J.  Math.  Anal. ,  19, 
172-203. 

[8]  Babuska,  I.,  Guo,  B.  Q.  [19881;  The  h-p  Version  of  the  Finite  Element 

Method  for  Domains  with  Curved  Boundaries,  SIAM  J.  Numer.  Anal.  25, 
837-861. 

[9]  Babuska,  I.,  Guo,  B.  Q.  [1989]:  The  h-p  Version  of  the  Finite  Element 

Method  for  Problems  with  Nonhomogeneous  Essential  Boundary  Conditions, 
Comput.  Methods  Appl.  Mech.  Engrg.  74,  1-28. 

[10]  Babuska,  I.,  Guo,  B.  Q.  [1987]:  The  Theory  and  Practice  of  the  h-p 

Version  of  the  Finite  Element  Method,  Computer  Method  in  Partial 
Differential  Equations  -  VI,  Ed.,  R.  Vichnevesky  &  R.  S.  Stepleman, 
241-247. 

[11]  Babuska,  I.,  Guo,  B.  Q.  [1989]:  Regularity  of  the  Solutions  of  Elliptic 

Problem  with  Piecewise  Analytic  Data,  Part  II.  The  Trace  Spaces  and 
Application  to  the  Boundary  Value  Problems  with  Nonhomogeneous  Boundary 
Conditions,  SIAM  J.  Math.  Anal.,  20,  763-781. 

[12]  Babuska,  I.,  Guo,  B.  Q. ,  Osborn,  J.  E.  [1989]:  Regularity  and  Numerical 
Solution  of  Eigenvalue  Problems  with  Piecewise  Analytic  Data.  SIAM  J. 
Numer.  Anal.  26,  1534-1560. 

[13]  Babuska,  I.,  Guo,  B.  Q. ,  Suri,  M.  [1989];  Implementation  of 
Nonhomogeneous  Dirichlet  Boundary  Conditions  in  the  p-Version  of  the 
Finite  Element  Method,  Impact  Comp.  Sci,  Engrg.  1,  36-63. 

[14]  Babuska,  I.,  Plank,  L. ,  Rodriguez,  R.  [1992]:  Quality  Assessment  of  the 


40 


A-posteriori  Error  Estimation  in  Finite  Elements.  Finite  Elements  in 
Analysis  and  Design.  To  appear. 

[15]  Babuska,  I.,  Rank,  E.  [1987]:  An  Expert  System  for  Optimal  Mesh  Design 
in  the  h-p  Version  of  the  Finite  Element  Method,  Internal .  J.  Numer. 
Methods  Engrg.  24,  2087-2106. 

[16]  Babuska,  I.  and  Rheinboldt,  W.  C.  [1978]:  A  Posteriori  Error  Estimates 
for  Adaptive  Finite  Element  Computations,  SIAM  J.  Numer.  Anal.,  15, 
736-754. 

[17]  Babuska,  I.,  Rheinboldt,  W.  C.  [1979]:  Analysis  of  Optimal  Finite 
Element  Meshes  in  R^ ,  Math.  Comp.  33,  435-463. 

[18]  Babuska,  I.,  Rheinboldt,  W.  C.  [1978]:  A  posteriori  Error  Estimators  in 

the  Finite  Element  Method,  Internal.  J.  Numer.  Methods  Engrg.  12, 
1597,1615. 

[19]  Babuska,  I.,  Rheinboldt,  W.  C.  [1980]:  Reliable  Error  Estimation  and 

Mesh  Adaptation  for  the  Finite  Element  Method,  in  Computational  Methods 
in  Nonlinear  Mechanics  (J.  T.  Oden,  ed.).  North  Holland  67-108. 

[20]  Babuska,  I.,  Suri,  M.  [1990]:  The  p-  and  h-p  Version  of  the  Finite 

Element  Method.  An  Overview.  Comput.  Methods  Appl.  Mech.  Engrg.  80, 
5-26. 

[21]  Babuska,  I.,  Suri,  M.  [1987].  The  h-p  Version  of  the  Finite  Element 

Method  with  Quasi-uniform  Meshes,  RAIRO,  Math.  Modelling  and  Numer. 

Anal.  21,  199-238. 

[22]  Babuska,  I.,  Szabo,  B.  A.,  Katz,  I.  N.  [1981]:  The  p-Version  of  the 
Finite  Element  Methods,  SIAM  J.  Numer.  Anal.,  18,  515-545. 

[23]  Babuska,  I.,  Vogelius,  M.  [1984]:  Feedback  and  Adaptive  Finite  Element 
Solution  of  One-dimensional  Boundary  Value  Problems,  Numer.  Math,  44, 
75-102. 

[24]  Babuska,  I.,  Yu,  D.  [1987]:  Asymptotically  Exact  a  Posteriori  Error 
Estimator  for  Biquadratic  Elements,  Finite  Elements  in  Analysis  and 
Design  3,  199-238. 

[25]  Bank,  R.  E. ,  Weiser,  A.  [1985],  Some  a-posteriori  Error  Estimators  for 
Elliptic  Partial  Differential  Equations,  Math.  Comp.  44,  283-301, 

[26]  Bank,  R.  E.  [1990]:  PLTMG:  A  Software  Package  for  Solving  Elliptic 
Partial  Differential  Equations.  User’s  Guide  6.0  SIAM  Philadelphia, 

PA. 

[27]  DaUge,  M.  [1988]:  Elliptic  Boundary  Value  Problems  on  Corner  Domains, 
Lecture  Notes  in  Math.  1341,  Springer,  New  York. 

[28]  Dorr,  M.  R.  [1984]:  The  Approximation  Theory  for  the  p-Version  of  the 
Finite  Element  Method,  SIMi  J.  Numer.  Anal.,  21,  1180-1207. 

[29]  Dorr,  M.  R.  [1986]:  The  Approximation  of  the  Solutions  of  Elliptic 
Boundary-value  Problems  Via  the  p-Version  of  the  Finite  Element  Method, 
SIAM  J.  Numer.  Anal.,  23,  58-77. 


41 


4 


[30]  Duran,  R. ,  Muschetti,  M.  A.,  Rodriguez,  R.  [1991]:  On  the  Asymptotic 
Exactness  of  Error  Estimators  for  Linear  Triangular  Elements,  Numer 
Math.  50,  107-127. 

[31]  Duran,  R. ,  Muschetti,  M.  A.,  Rodriguez,  R.  [1992]:  Asymptotically  Exact 
Error  Estimators  for  Rectangular  Finite  Elements,  SIAM  J.  Numer.  Anal., 
28,78-88. 

[32]  Duran,  R. ,  Rodriguez,  R.  [1991):  On  the  Asymptotic  Exactness  of 
Bank-Weiser  Estimators,  Numer.  Math.  To  appear. 

[33]  Eriksson,  K. ,  Johnson,  C.  [1988]:  Adaptive  Finite  Element  Method  for 
Linear  Elliptic  Problems,  Math.  comp.  50,  361-368. 

[34]  Ewing,  R.  E.  [1990]:  A  posteriori  Error  Estimation,  Comput.  Methods 
Appl.  Mech.  Engrg. ,  82,  59-72. 

[35]  Grisvard,  P.  [1985]:  Elliptic  Problems  in  Nonsmooth  Domains,  Pitman, 
Boston. 

[36]  Gui,  W. ,  Babuska,  I.  [1986]:  The  h-p  Versions  of  the  Finite  Element 
Method  in  One  Dimension,  Part  1:  The  Error  Analysis  of  the  p-Version; 
Part  2:  The  Error  Analysis  of  the  h  and  h-p  Versions,  Part  3:  The 
Adaptive  h-p  Version,  Numer.  Math. , A3,  577-612,  613-657,  659-683. 

[37]  Guo,  B.  Q.  [1988]:  The  h-p  Version  of  the  Finite  Element  Method  for 
Elliptic  Equation  of  Order  2m,  Numer.  Math.  53,  199-224. 

[38]  Guo,  B.  Q. .Babuska,  I.  [1986[:  The  h-p  Version  of  the  Finite  Element 
Method,  Part  1:  The  Basic  Results,  Part  2:  General  results 

and  Applications,  Comput.  Mech.  1,  21-41,  203-230. 

[39]  Hugger,  J.  [1992]:  The  Theory  of  Density  Representation  of  Finite 
Element  Meshes.  Examples  of  Density  Operators  with  Quadrilateral 
Elements  in  the  Mapped  Domain.  To  appear. 

[40]  Hugger,  J.  [1992];  Recovery  and  Few  Parameter  Representations  of  the 
Optimal  Mesh  Density  Function  for  Nearly  Optimal  Finite  Element  Meshes. 
To  appear. 

[41]  Kelly,  D.  W.  [984]:  The  Self  Equilibration  of  Residuals  and 
Complementary  Error  Estimates  in  the  Finite  Element  Method,  Internet.  J. 
Numer.  Methods  Engrg. ,  20,  1491-1506. 

[42]  Kondratev.  V.  A.  [1967]:  Boundary  Value  Problems  for  Elliptic  Equations 
in  Domains  with  Conical  or  Angular  Points,  Trans.  Moscow  Math.  Soc.  16, 
227-313. 

[43]  Maday,  Y. ,  Patera,  A.  T.  [1989]:  Spectral  Element  Methods  for  the 
Incompressible  Navier-Stokes  Equation  in  State  of  the  Art  Surveys  on 
Computational  Mechanics,  A.  K.  Noor,  J.  T.  Oden  (eds. ]  Amer.  Soc.  Mech. 
Engr. ,  New  York,  71-143. 

[44]  Mestenyi,  C. ,  Szymczak,  W.  [19811;  FEARS  user’s  Manual  for  UNIVAC  1100. 
University  of  Maryland,  Institute  for  Physical  Science  and  Technology, 


42 


4 


Technical  Note  BN-991. 

[45]  MSC/PROBE  [1990],  User’s  manual.  MacNeal  Schwendler  Co.  Los  Angeles,  CA. 

[46]  Noor,  A.  K.  ,  Babuska,  I.  {19871;  Quality  Assessment  and  Control  of 
Finite  Element  Solutions,  Finite  Elements  in  Analysis  and  Design  3, 

1-26. 

[47]  Oden,  J.  T.  [1988]:  Adaptive  Finite  Element  Method  for  Problems  in 
Solid  and  Fluid  Mechanics,  in  Finite  Elements,  Theory  and  Application, 

D.  L.  Dwyer,  M.  Y.  Hussani,  R.  B.  Voigt  (eds.).  Springer,  268-291. 

[48]  Oden,  J.  T.  ,  Demkowicz,  L. ,  Rachowicz,  W. ,  Westermann,  T.  A.  [1989], 
Towards  a  Universal  h-p  Finite  Element  Strategy:  Part  2,  A  posteriori 
Error  Estimation,  Compat.  Methods  Appl.  Mech.  Engrg.  77,  113-180. 

[49]  Oden,  J.  T. ,  Rachowicz.  W,  Demkowicz,  L.  [1989]:  Toward  a  Universal  h-p 
Adaptive  Element  Strategy,  HI.  Design  of  h-p  meshes,  Comput.  Methods 
Appl.  Mech.  Engrg.  H,  181-1212. 

[50]  PHLEX  User’s  Manual  [1990],  Computational  Mechanics  Co.,  Austin,  TX. 

[51]  Rivara,  M.  C.  [1984];  EXPDES  User’s  Manual,  Catholic  University, 

Leuven,  Belgium. 

[52]  Shephard,  M.  S. ,  Finnigan,  P.  M.  [1989]:  Toward  Aautomatic  Model 
Ggeneration  IN  State-of-the-Art  on  Computational  Mechanics,  (A. 

K.  Noor  and  J.  T.  Oden,  eds.),  Amer.  Soc.  Mech,  Engr. ,  New  York. 

[53]  Shephard,  M.  S. ,  Niu,  Q. ,  Bachmann,  P.  L.  [1989];  Some  Results  Using 
Stress  Projectors  for  Error  Indication  and  Estimation,  in  Adaptive 
Methods  for  Partial  Differential  Equations,  J.  E.  Flaherty,  P.  J. 

Paslow,  M.  S.  Shephard,  J.  D.  Vasilakes,  eds. ,  SIAM  Philadelphia, 

83-99. 

[54]  Szabo,  B.  A.  [1986]:  Estimation  and  Control  of  Error  Based  on 
P-convergence ,  in  Accuracy  Estimates  and  Adaptive  Refinements  in  Finite 
Element  Computations,  (I.  Babuska,  0.  C.  Zienkiewicz,  J.  Gago,  E.  R.  de 
Olivera  eds.),  J.  Wiley  &  Sons,  New  York,  61-78. 

[55]  Szabo,  B.  A.  [1990]:  The  p  and  h-p  Versions  of  the  Finite  Element 
Methods  in  Solid  Mechanics,  Comput.  Methods  Appl.  Mech.  Engrg.  80, 
185-195. 

[56]  Szabo,  B.  A.,  Babuska,  I.  [1991]:  Finite  Element  Analysis,  J.  Wiley  & 
Sons,  New  York. 

[57]  Zienkiewicz,  0.  C. ,  Zhu,  J.  Z.  [1987]:  A  Simple  Error  Estimator  and 
Adaptive  Procedure  for  Practical  Engineering  Analysis,  Internal .  J. 

Numer.  Methods  Engrg.  24,  337-357. 

[58]  Zienkiewicz,  0.  C. ,  Zhu,  J.  Z.  [1990];  The  three  R’ s  of  engineering 
analysis  and  error  estimation  and  adaptivity,  Comput.  Methods.  Appl. 

Mech.  Engrg.  82,  95-113. 

[59]  Zienkiewicz,  0.  C.  Zhu.  J.  Z. ,  Gong,  N.  G.  [1989]:  Effective  and 
Practical  h-p  Version  Adaptive  Analysis  Procedures  for  Finite  Element 


43 


Method,  Internal.  J.  Numer.  Methods  Engrg.  28.  879-891. 

[60]  Zienkiewicz,  0.  C.  Zhu.  J.  Z.  [1992]:  Superconvergent  Derivative 
Recovery  Techniques  and  a  Posteriori  Error  Estimation  in  the  Finite 
Element  Method,  Part  1:  A  general  superconvergent  recovery  technique. 
Report  of  INME,  University  College  of  Swansea,  CR/671/91,  also  Internal. 
J.  Numer.  Methods  Engrg.  To  appear. 

[61]  Zhu.  J.  Z. ,  Zienkiewicz,  0.  C.  [1990]:  Superconvergence  Recovery 
Technique  and  a  posteriori  Error  Estimators,  Internal  J.  Numer.  Methods 
Engrg.  30,  1321-1339. 


44 


The  Laboratory  for  Numerical  Analysis  is  an  integral  part  of  the  Institute  for  Physical 
Science  and  Technology  of  the  University  of  Maryland,  under  the  general  administration  of  the 
Director,  Institute  for  Physical  Science  and  Technology.  It  has  the  following  goals: 

•  To  conduct  research  in  the  mathematical  theory  and  computational  implementation  of 
numerical  analysis  and  related  topics,  with  emphasis  on  the  numerical  treatment  of 
linear  and  nonlinear  differential  equations  and  problems  in  linear  and  nonlinear  algebra. 

•  To  help  bridge  gaps  between  computational  directions  in  engineering,  physics,  etc.,  and 
those  in  the  mathematical  community. 

•  To  provide  a  limited  consulting  service  in  all  areas  of  numerical  mathematics  to  the 
University  as  a  whole,  and  also  to  government  agencies  and  industries  in  the  State  of 
Maryland  and  the  Washington  Metropolitan  area. 

•  To  assist  with  the  education  of  numerical  analysts,  especially  at  the  postdoctoral  level, 
in  conjunction  with  the  Interdisciplinary  Applied  Mathematics  Program  and  the 
programs  of  the  Mathematics  and  Computer  Science  Departments.  This  includes  active 
collaboration  with  government  agencies  such  as  the  National  Institute  of  Standards  and 
Technology. 

•  To  be  an  international  center  of  study  and  research  for  foreign  students  in  numerical 
mathematics  who  are  supported  by  foreign  governments  or  exchange  agencies 
(Fulbright,  etc.). 

Further  information  may  be  obtained  from  Professor  I.  Babuska,Chairman,  Laboratory  for 
Numerical  Analysis,  Institute  for  Physical  Science  and  Technology,  University  of  Maryland,  College 
Park,  Maryland  20742-2431. 


